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Abstract. We propose an accurate tight-binding parametrization for the band structure of 
M 0 S 2 monolayers near the main energy gap. We introduce a generic and straightforward 
derivation for the band energies equations that could be employed for other monolayer 
dichalcogenides. A parametrization that includes spin-orbit coupling is also provided. The 
proposed set of model parameters reproduce both the correct orbital compositions and location 
of valence and conductance band in comparison with ab initio calculations. The model gives 
a suitable starting point for realistic large-scale atomistic electronic transport calculations. 


1. Introduction 

The synthesis of graphene in 2004 null, the first single-atom thick material, has boosted 
the research in atomically thin two-dimensional ( 2 D) materials. The ability to manipulate 
isolated single atomic layers and reassemble them to form heterostructures layer-by-layer 
in a precise sequence, opens enormous possibilities for applications El a 12 a. Along 
this approach, semiconducting dichalcogenides are promising compounds since they can 
be easily exfoliated and present a suitable small gap both in bulk and as a single layer. 
In this category of 2 D dichalcogenides systems, monolayer molybdenum disulfide (M0S2) 
has recently gained attention for combining an electron mobility comparable to that of 
graphene devices with a finite energy gap Q. Unlike its bulk form, which is an indirect 
gap semiconductor, monolayer M0S2 has a direct gap Eia, making it very interesting for 
optoelectronics. Another interesting feature is that the electronic properties appear to be 
highly sensitive to external pressure strain ifTOl im . and temperature Ea, which affect 
the gap and, under certain conditions, can also induce a insulator/metal transition. In addition, 
the lack of lattice inversion symmetry together with spin-orbit coupling (SOC) leads to 
coupled spin and valley physics in monolayers of M0S2 and other group-VI dichalcogenides 
M M, making it possible to control spin and valley in these materials El El. Due 
to their peculiar band structure, a variety of nanoelectronics applications Elia including 
valleytronics, spintronics, optoelectronics, and room temperature transistor devices jT) have 
been suggested for monolayers of M0S2. 

In light of the growing interest in this material, an accurate and yet reasonably simple 
model describing the band structure and electronic properties of M0S2 is highly desirable. So 
far the electronic properties of single-layer and few-layer dichalcogenides have been mainly 
investigated by means of ab initio calculations, based on Density Functional Theory (DFT) 
Einni. Such methods provide valuable information about electronic properties of pristine 
dichalcogenide crystals, but are computationally prohibitive to treat disordered systems with 
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a large number of atoms. To address the latter, one needs to resort to a simple effective model, 
such as the kp Hamiltonian or the tight-biding approximation. In this paper we choose the 
latter route, which provides a more accurate description for the entire band structure than 
the kp method. Moreover, the tight-binding model applied to a single-layer M 0 S 2 as well 
as to similar transition metal dichalcogenides, constitutes a key tool for further studies of 
the low-energy electronic transport properties of these materials, such as the description of 
the conductivity in diffusive samples, as well as for evaluating the conductance of ballistic 
samples as a function of carrier concentration. 

In recent years, a variety of tight-binding models have been proposed for M 0 S 2 
monolayers ESI El [m. Unfortunately, they are neither practical nor sufficiently accurate for 
transport calculations. For that purpose, one needs a tight-binding model with a manageable 
number of parameters and interactions that accurately reproduces the ab initio electronic 
properties of the conduction band (CB) near its maximum points and the the valence band 
(VB) near its minimal points. Before we present our results, let us now briefly review the 
main features of the tight-binding models for dichalcogenides found so far in the literature. 

An “all orbital model” was put forward by Zahid and collaborators ifTbll . The model 
includes non-orthogonal sp^d^ orbitals, considers only nearest-neighbour hopping matrix 
elements, and includes spin orbit coupling. The model has 96 fitting parameters. The 
optimization of the Slater-Koster energies ED and overlap integrals used in the model are 
obtained by a fit to the DFT target band structure. The model shows good agreement with 
band structure calculations using the HSE06 functional EOl 1^ . but its computational cost 
and complexity make it impractical for studying disorder and electronic transport at large 
scales. 

In contrast, Liu and collaborators ED proposed a three-orbital tight-binding model. 
The authors consider a superposition of orbitals dp.,dxy, and dyi_y 2 as orthogonal basis, 
targeting the main orbital composition around the K point, which corresponds to the direct 
gap. Thus, the agreement between their first nearest-neighbour tight-binding model and 
the DFT predictions using both the local density approximation (LDA) ll22l [23 and the 
generalized-gradient approximation (GGA) ll2^ is limited to features in the vicinity of the 
K point, missing the local band minimal at the Q point. By including up to the third-nearest 
neighbour hopping involving Mo-Mo terms, the agreement with the DFT-GGA band structure 
improves substantially. This is achieved at the expense of increasing the complexity of the 
model, as the number of fitting parameters goes from 8 to 19. In transport calculations, the 
inclusion of higher neighbour hopping terms implies in an increase in the size of the unit cell. 
Hence, trading a larger number of bands with nearest-neighbour hopping for a simpler model 
with longer range hopping is not necessarily advantageous. Moreover, the orbital composition 
in Ref. EH is, by construction, restricted to Mo orbitals which limits the analysis of disorder 
effects. It is also worth noting that Ref. Ell fails to reproduce the orbital composition and 
energy spectrum around the F point, which plays a significant role in transport for hole-doped 
monolayers. 

A seven-orbitals tight-binding parametrization has been introduced by Rostami and 
collaborators ESl . The model considers a non-orthogonal basis and neglects the s and p^ 
orbitals of the S atoms and the s, dy^, and dxz orbitals of the Mo atom by invoking arguments 
based on crystal symmetry and the range of energies of interest. The model reproduces the 
main features around the K point, but two unrealistic flat bands appear in the gap region. We 
attribute this undesired feature to the fact that the p^ orbital of the S atoms is not actually 
decoupled from other orbitals and can not be neglected. On the contrary, the p^ orbital from S 
atoms plays a pivotal role in the transition from a direct to an indirect gap, when passing from 
a monolayer to a multilayer system. The basis set introduced in Ref. ESl does not distinguish 
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between p[, and S orbitals at the top and bottom planes of the S-Mo-S layers. Therefore, 
it can not correctly capture the symmetry under inversion of the z-axis. 

As pointed out by Cappelluti and collaborators ESll . the linear combination of p* and 


p^ orbitals is necessary to produce z-symmetric and z-antisymmetric states. For this reason. 
Refs, ll^ l26l propose a minimal model of 11 orbitals. This is also our choice. This model 
considers an orthogonal basis composed of all the Ad Mo orbitals and the 3p S orbitals, 
forming real symmetric (even) and antisymmetric (odd) combinations of Px,y,z orbitals 
lIZTl l28l . The tight-binding parameters found in Refs. Il25l l26ll yield two bands that look 
very similar to the conduction and the valence bands obtained by standard DFT calculations. 
However, in our treatment, by using analytical expressions for the valence and conductance 
bands at high-symmetry k-points, we observe that the tight-binding orbital compositions of 
Refs. 12511261 have actually no relation with those calculated using DFT. Hence, a new tight- 
binding parametrization, reproducing both energies and orbital composition is badly needed. 
This is the main goal of this paper. We rederive the tight-binding equations of Ref. 1251 
in a more direct and transparent way, allowing us to more carefully consider the orbital 
composition in our parameter optimization procedure. 

The paper is organized as follows. In Sec.|^we describe the atomic structure of a M 0 S 2 
monolayer and discuss the DFT-HSE06 results for the the band structure that will be the 
reference for our tight-binding model. In Sec. we present the model. In Sec. we analyze 
the band equations for a few high-symmetry k-point, allowing us to obtain simple analytical 
expressions for the bands. These are used to hnd the best set of tight-binding parameters 
that ht the DFT band structure. In Sec. we present the optimized parameters and the 
corresponding band structure. In Sec.|^we consider a simplihed model with a reduced number 
of parameters. In Sec. we add spin-orbit interaction to the full model. Finally, in Sec. we 
draw our conclusions. 

The main text is supplemented by a number of appendices containing technical aspects of 
the calculations. In [Appendix A we present in detail all the elements required to construct the 
tight-binding band equations. In [Appendix Bj we show how to implement the band equations 
in the unsymmetrized and symmetrized ones. Appendix C analyzes the band structure at 
symmetry points used in the optimization. Finally, in Appendix D we present a comparison 
between our 11-band tight-binding formulation and that of Refs. l25]|26l. 


2. Crystal structure and ab-initio electronic structure 


Molybdenum disulhde is a layered transition metal dichalcogenide semiconductor. The 
layered structure is formed by a honeycomb arrangement of Mo and S atoms stacked together 
and forming S-Mo-S sandwiches coordinated in a triangular prismatic fashion. The S-Mo-S 
layers are bonded together by weak van der Waals forces. 

The single-layer M 0 S 2 lattice structure is shown in Fig.[2 top and lateral views. It is a 
2D rhombic lattice with a three-atom basis (one Mo and two S). The two Bravais primitive 
lattice vectors are 


and 


Ri = (a,0,0) 


R2 = 



( 1 ) 

( 2 ) 


where a = 3.16 A is the lattice constant. The S atoms are located in planes 1.56 A above 
and below the Mo plane. This yields a distance between neighboring Mo and S atoms of 
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d = 2.40 A. The angle between the Mo-S bond and the Mo plane is 9b — 40.6". These values 
are obtained by the DFT calculation discussed below and are consistent with previous DFT 
calculations osiiioiiia and with experimental values Esma. 


z 




Figure 1: Left panel: Top view of the M 0 S 2 lattice structure. Dark (light) circles represent 
Mo (S) atoms. Notice that in this view two S atoms sit on top of each other. The unit cell is 
shown in the highlighted hexagon. The lattice constant in the Mo plane is a. The two Bravais 
lattice vectors (Rj and R 2 ) are indicated. Six other auxiliary vectors that connect a Mo atom 
with its nearest S atoms, 5i±, 82 ±, and 53±, are indicated. Right panel: Tridimensional view 
of the first neighbors of a Mo atom. The reference trigonal prism coordination unit and other 
useful quantities are also shown. 


For the purpose of building the tight-binding model, we will follow the notation 
introduced Fig. We denote by “f” (or “+”) and by “Z?” (or the S atoms at the top 
and bottom layers, respectively. The distance between the two S layers is dcosOg = ajsffi. 
The nearest-neighbour vectors, connecting Mo and S atoms, are given by 


5i± = £/ (OjCos 0B,±sin0B), 

82 ± =d I-^cos0B,-icos0B,±sin0B I , 


53± =d I -i-^COS0B,-^COS0B,±sin0B 


(3) 

(4) 

(5) 


The M 0 S 2 Brillouin zone is hexagonal. The most important symmetry points and symmetry 
lines are indicated in Fig. namely, F = (0,0), K = , and M = ■ The 

reciprocal lattice basis vectors are 


Ki 


471 f V3 
s/3a y 2 



(6) 


and 

47r 

K2 = ^ (0,1,0) (7) 

\'3a 

Table [ 1 ] summarizes the experimental and theoretical values of the band gap of M 0 S 2 . 
Early photoluminescence experiments [31291 had inferred a direct band gap of about 1.9 eV 
for MoS 2 . More recently, it has been suggested that this value is actually the result of excitonic 
states and hence corresponds to the optical gap rather than the actual direct gap between 
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Figure 2: Brillouin zone for the M 0 S 2 lattice. Ki and K 2 are the reciprocal lattice basis 
vectors, and F, K, K', and M are the high-symmetry points considered in this study. 

Table 1: Summary of experimental and theoretical values of the band gap of M 0 S 2 . 


Method 

Value (eV) 

Note 

Photoluminescence 1811291 

1.8-1.9 

Optical band gap 

Scanning Tunneling Spectrocopy (STS) ETI 

2.15 


DET-LDA 

1.81 


DET-PBE Eli 

1.68 


DET-optB88-vdW ES 

1.67 


GW l30l 

2.84 

GiWq approximation 

DET-HSE06 

2.23 

This work 


the single-particle VB and CB ll^ . Scanning Tunneling Spectroscopy (STS) measurements 
revealed that the band gap of M 0 S 2 is 2.15 eV ED- Given the optical gap of about 1.9 
eV, the latter value is quite consistent with both theoretical and experimental values of the 
exciton binding energy, which fall in the range 0.28-0.33 eV according to theory ll^ 
and are either 0.44 eV Ga or 0.22 eV 1311 as deduced from experiments. Traditional DFT 
functionals based on the local density approximation (LDA) and on the generalized gradient 
approximation (GGA), not surprisingly, underestimate this band gap llTSl lJ7l . while the 
more advanced GW approach tends to overestimate it OOl . The HSE06 functional 1201IZD, 
on the other hand, provides so far the best agreement 1381 with the STS result for this gap 

ED- 

In this work, we have therefore chosen the DFT-HSE06 band structure as reference for 
our htting procedures. Our DET-based electronic band structure calculations are carried with 
the HSE06 functional using the supercell method with a plane-wave basis set (cutoff energy 
of 500 eV) and the projector-augmented wave (PAW) technique l39l l40l . as implemented 
in the Vienna ab-initio Simulation Package (VASP) EUia. We use a supercell consisting 
of a M 0 S 2 layer with an experimental lattice parameter value of 3.16 A at its center and 
a vacuum of 15 A to minimize the interaction between normal periodical images. The 
structure is optimized using the GGA approximation with the Perdew-Burke-Ernzerhof (PBE) 
parameterization l24l . The Brillouin zone is sampled byal8xl8xl mesh. In calculations 
including spin-orbit-coupling, we sample the Brillouin Zone with a 9 x 9 x 1 k-point mesh to 
reduce computational cost. The electronic band structure along the Y-K-M-Y directions is 
calculated with 149 k-points and then projected onto every orbital of each atom to resolve the 
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symmetry characters of the cotTesponding wave-functions. The resulting band shvicture for a 
M 0 S 2 monolayer is shown in Fig. 



Figure 3: The DFT-HSE06 band structure of M 0 S 2 near the gap region. See text for details. 
Let us summarize the main features near the high-symmetry points of the Brillouin zone 

MM: 

• K point - The calculated (DFT-HSE06) band gap Eg = 2.23 eV is located at the K point. 
This result is good agreement with experimental value of 2.15 eV Ell- Electron-hole 
symmetry is clearly absent: While the effective mass in the CB is nearly isotropic, in 
the VB band it is characterised by a strong trigonal warping. Spin splitting is present in 
both CB and VB, but the splitting of the VB states is much larger. The VB spin splitting 
at the K point is experimentally found to be around 145 meV. Here we only consider 
spin-orbit coupling up to hrst order in the coupling and hence disregard the spin splitting 
in the CB at the K point. Higher order SOC effects have been studied in Refs. 14411451 . 

• Q point - This point, signaling a local minimum in the CB band along the straight line 
connecting F and K points, has recently received increased attention due to its relevance 
for transport properties, since the energy minimum Eq is close to the bottom of the 
CB ESlIll. From our DFT-HSE06 calculations, before including SOC, we estimate 
this energy difference to be AE « 0.3 eV. It is noteworthy that the CB at the Q point 
moves down in energy in multilayer systems. As discussed in Ref. Eh), phonon-limited 
mobility depends quite sensitively on this energy separation. At the Q point, the CB 
is characterized by a spin splitting of 91 meV and the effective mass has an ellipsoidal 
shape US). The Q point is located close to the mid point between F and K points. 
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• r point - This point lies close to the top of the valence band. According to our DFT- 
HSE06 calculations, before including SOC, its energy difference to the K point is very 
small, namely, AE « 0.15 eV. Hence, in hole-doped samples states both the K and F 
points will contribute to the electronic transport. 

The orbital composition is of fundamental importance for building of any tight-binding 
model. As already found in literature ll25l . dxy,d^ 2 _y 2 ,dp, and pxj are the most important 
orbitals to describe the valence and the conduction bands. It is worth to stress that dy^ and d^^ 
have a dominant contribution at the F point, and gives an important contribution to both F 
and Q points. Thus, for a comprehensive description of the CB and VB along the Brillouin 
zone, one needs to consider all these orbitals. Tables]^ andshow the relative contribution 
from each orbital at the points F and K, as provided by our DFT-HSE06 calculations. These 
results serve not only to justify the choice of the relevant orbitals of the atomistic model, but 
also help in finding the right constraints to optimize the tight-binding parameters. 

Table 2: Density functional theory (DET-HSE06) orbital composition at the F point for 
different bands. Absence of an entry indicates zero contribution. 


band 

number 

band 

energy (eV) 

Py 

Pz 

Px 

^xy 

dyz 

dz2 

dxz 

dy2 

6 

-7.571 


0.252 




0.199 



7 

-4.105 

0.115 


0.127 

0.309 




0.279 

8 

-4.105 

0.127 


0.115 

0.279 




0.309 

9 

-3.303 


0.558 







10 

-2.753 

0.009 


0.313 


0.015 


0.502 


11 

-2.753 

0.313 


0.009 


0.502 


0.015 


12 VB 

-1.262 


0.141 




0.596 



13 CB 

2.457 

0.302 


0.045 


0.327 


0.049 


14 

2.457 

0.045 


0.303 


0.048 


0.327 


15 

2.678 

0.083 


0.390 

0.249 




0.053 

16 

2.678 

0.390 


0.083 

0.053 




0.249 


3. Model 

Let r, denote the Mo atom location in the /th unit cell. Eollowing Cappelluti and collaborators 
Il25ll2^ . we consider a tight-binding model with five d orbitals in the Mo atom, namely, 

|r;;t/o) = M3 z2-f). = \dx2^y2), |r;;c/2) = \dxy), |r/;t/3) = \rf,d4) = \dyf), 

and six p orbitals for the S atoms, three for the top t (-f) and three for the bottom b (—) layers, 
\ri + 5i±;pi) = \p‘/), |ri + 5i±;p2) = IK’*), |d + 5i±;/73) = |K*)- 

Starting with this basis we can define on-site energies and hopping amplitudes and write down 
a tight-binding Hamiltonian. Hereafter we assume that this basis is orthogonal. 

The tight-binding Hamiltonian Jff contains Mo-S and S-S nearest-neighbor hopping 
terms (in the same unit cell), as well as Mo-Mo and S-S next-to-nearest-neighbor ones (in 
adjacent cells). Each Mo has six S nearest neighbors, while the next-to-nearest neighbor 
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Table 3: Density functional theory (DFT-HSE06) orbital composition at the K point for 
different bands. Absence of an entry indicates zero contribution. 


band 

number 

band 

energy (eV) 

Py 

Pz 

Px 

dxy 

dy^ 

d^2 

dxz 

d^i 

6 

-7.259 

0.155 


0.155 


0.184 


0.184 


7 

-6.427 

0.178 


0.178 



0.135 



8 

-5.742 

0.231 


0.231 

0.034 




0.034 

9 

-5.244 

0.034 

0.359 

0.034 

0.105 




0.105 

10 

-4.466 

0.230 

0.129 

0.230 






11 

-3.734 

0.416 



0.140 


0.140 



12 VB 

-1.111 

0.065 


0.065 

0.345 




0.345 

13 CB 

1.120 

0.034 


0.034 



0.753 



14 

2.718 

0.068 


0.068 


0.248 


0.248 


15 

3.284 

0.016 

0.153 

0.016 

0.327 




0.327 

16 

4.899 


0.188 



0.303 


0.303 



hoppings connect 6 atoms of the same kind, see Fig.[2 Overall, there is a total of 25 hopping 
matrix elements inside the unit cell and between the unit cell and the adjacent cells. 

The hopping amplitudes are written in terms of Slater-Koster (SK) parameters QSl. We 
incorporate the x and z reflection symmetries in the construction of the basis, when applicable. 


to reduce the number of terms. We refer to Appendix A for a detailed description of the tight- 
binding Hamiltonian and the transfer integrals. There, we also provide expressions for the 
hopping amplitudes in terms of the SK integrals Vp^a,Vpdjc,Vdda,VddS,Vddn,Vppa, and Vpp^- 
This allows for a significant reduction in the number of fitting parameters of the model. 

To And the energy bands we solve the eigenvector equation that, in the Bloch momentum 
representation, reads 


Jff’lk) =£a(k)|k) 


( 8 ) 


|k)=E.'"--- 


where the eigenstates |k) are expressed in terms of the three-atom basis, namely, 

4 

.M=0 
3 

E {Pk.n\ri + Si--,pp,) + Tk^pa\ri + Si+-,pp)) 

M=1 

For the purpose of implementing the eigenvalue equation we project the vector ,^|k) 
onto the three-atom basis and write 

Jt \ / a\ /a 

{py h^ + U S 1 I T 1 =£ I T I , (10) 


(9) 


(r 


b\f 


h^ + U 


where have omitted, for the moment, the spin i ndices. Expli cit expressions for the block 
matrices h^°,h^,T',T'’,S,U, and V are given in [Appendix A[ The matrices S, U and V are 
real and symmetric. 

The secular equation ( [T0| ) is sufficient for a numerical evaluation of the band structure. 
Nonetheless, it is convenient to use symmetry arguments to reduce the the size of the matrices 
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to be diagonalized, allowing us to express analytically the gap and other features of the band 
structure at k points of interests. 

Let us introduce the symmetric and anti-symmetric components 


and 


0k,V — ("^k.v + ^k,v) 


^k,v — ~ Pk.v)■, 


( 11 ) 


( 12 ) 


(13) 


that allow us to write the Hamiltonian in the matrix form 

h^° \ / a \ (a 

h^P 0 0 = £ 0 

jO-i Q J \ (j) J \ (j) 

The eigenvalue problem can be further simplified by rearranging rows and columns through 
the transformation \j/ \(r, where 

V/^ = (ao,ai,a2,o:3,a4,0i,02,03,0i,02,03) (14) 

and 

if = (ao,ai,a2,0i,02,03,o!3,a4,0i,<)>2,03)- (15) 

Notice that the first six orbital basis functions are even (E) with respect to a z-axis 
inversion, while the last five are odd (O). Then, the problem is reduced to two decoupled 
eigenvalue/eigenvector problems, namely 

0 

0 ho 


\j/ = E\j/. 


(16) 


We refer to 


Appendix B|for explicit expressions of the matrix elements of and H^. 


4. Optimization of model parameters 


Our tight-binding model Hamiltonian has Np = 12 fitting parameters, namely, five on¬ 
site orbital energies {DQ,Di,D 2 ,Dp, and D^) and seven SK parameters related to hopping 
(ypdn,ypda,yppa,Vppn,Vdda,Vddx, and These parameters are optimized to reproduce 

the main characteristics of the low-energy bands we obtained from DFT-HSE06 calculations. 

Our main goal is to reproduce the energies, orbital composition, and effective masses 
of the conduction and valance bands at the 2, and T points. For that purpose we choose 
a number of representative A:-points, shown in Fig. and collect the corresponding band 
energies £’«(k), where n is the band index, to built the data set to be fitted. To better describe 
the conduction and valence energy bands, we give a larger weight to points (k,£’„(k)) near 
the main band gap. In addition, we take a larger concentration of points around K,Q, and F 
to reproduce the electron effective mass around these high symmetry-points. 

We find the optimal tight-binding parameters using the method of least squares. The data 
set is built from the band energies E^^, where j labels both the k-point and band index (see 
Fig.|^, with j = 1, • • • ,Nf. The corresponding j^-squared function is just a sum of weighted 
squared residuals, namely, 

2 


Nf 

5(P) = t 

j=i 


(P)-E 


DPT 


O’? 


( 17 ) 
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> 
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r Q K M r 

Figure 4; Reference DFT-HSE06 band structure with the constraint points indicated. Blue 
circles; analytical constraints. Orange circles: numerical constraints. Predominantly even 
(odd) bands with respect to z inversion are shown in blue (red). 
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4.0 

2.0 

0.0 

- 2.0 

-4.0 

- 6.0 


where (7, is the weight given to the y'th data set fitting point, £*(?) is the tight-binding energy 
corresponding to the yth data set fitting point. The vector array P of dimension Np contains 
the tight-binding parameters to be minimized. We minimize S with respect to P using the 
Powell method Hbl . that is an efficient method to hnd the minimum of a function of several 
variables without requiring the computation of its derivatives. 

Let us briefly describe the route we follow to approximate the low-energy band structure 
by progressively adding data points. 


(i) We compare the results obtained at the K and F points using the analytical expressions 
derived in Appendix C|to the DFT-HSE06 energy values and their orbital compositions. 


(ii) We consider k points in the vicinity of K and F. The weights a are adjusted to decrease 
the importance of these k points as the further away they are from the gap region. 


(iii) We consider additional data points to correctly reproduce other features of the CB and 
VB. In particular, we add k points at and around the M and Q points, to obtain a correct 
band energy behavior around the main gap over the entire Brillouin zone. 


The tight-binding model Hamiltonian decouples into “even” bands (associated to He) 
and “odd” bands (associated to Hq)- The identification of the parity of the bands at the F and 
K points allows us to follow all the bands over the entire Brillouin zone. Notice that around 
the main energy gap, bands are mostly even (except for the CB at the F point). Let us now 
explain how we match DFT-HSE06 band energies with the tight-binding eigenvalues at the K 
and F-points. 


4.1. K point 


At the K point, the matrices and in Eq. (16 1 can be written in block diagonal form, 
namely, we can break into three 2x2 diagonal blocks and into two 2x2 blocks and 
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one 1x1 block. The explicit expressions are presented in [Appendix C| 

The inspection of the orbital composition given by the DFT-HSE06 calculations (see 
Table allows us to establish a unique correspondence between pairs of bands and the 2x2 
blocks mentioned above. The correspondence is summarized in Table 

The identihcation of the highest £+ (lowest £_) eigenvalue of a given block with the 
highest (lowest) energy value of the band with corresponding orbital composition and parity 
severely constrains our model parameters. By applying this methodology to all diagonal 
blocks, we find analytical expressions for the band energies at the K point, presented in 
[Appendix Cj 

Table 4; Identification of the tight-binding 2x2 block structures and their orbital contributions 
at the K point with the band numbers and their corresponding DFT-HSE06 energies, given in 
Table 


block structure orbitals 

orbital composition 

(£+,£_) (eV) 

band numbers 

He : cxq — Or 

d^2,Px,Py 

(1.120,-6.427) 

(13CB, 7) 

He : ccli — Or 

d-x^ dxy tPxtPx 

(-1.111,-5,742) 

(12VB, 8) 

He : ctRi — 03 

dx^-y^^ dxy, Pz 

(3.284,-5.244) 

(15,9) 

Ho '■ (Xri — 0L 

dyz, dxz, Px,Py 

(2.457,-2.753) 

(14,6) 

Ho ■ 0Cl\ - O3 

dyz, dxz, Pz 

(4.899,-3.734) 

(16,11) 

Ho '■ 

Px,Py 

-4.466 

10 


4.2. r point 

At the r point we can also express the H matrix of Eq. ( fTbj ) in block block diagonal form, 
breaking it into hve 2x2 blocks and one 1x1 block. 

The 


Here we follow the procedure described in the previous subsection. Sec. 4.1 


differences are due to the distinct point-group symmetries of the K and E points. As a 
consequence, the orbital compositions of the tight-binding 2x2 blocks considered in this case 
are not the same as for the K point. This issue is discussed in [Appendix C[ where we also 
present the analytical derivation of eigenvalues and eigenstates at the E-point. 

Table [^presents the identihcation of the tight-binding symmetry split 2x2 blocks with 
their corresponding DFT-HSE06 bands. It is also noteworthy that, as presented in Table|^ the 
ad initio calculations show that several band energies coincide at the E point, namely, 7 and 
8 , 10 and 11, 13CB and 14, and 15 and 16. 


5. Eleven-band model: parameters and results 

In this Section we present the main results of our study, namely, the tight-binding 11-band 
parametrization and the corresponding band structure for M 0 S 2 . Table presents the best 
htting parameters we obtained using the the optimization procedure described in Sec.|^ 

Before discussing the results, it worth mentioning that the even-odd parity symmetry of 
our tight-binding model prevents a perfect match with ab initio calculations. For instance, 
DFT-HSE06 calculations indicate that the CB and VB are mainly “even”, but around the E 
point they gain a significant odd contribution. Despite this proviso, we show that the tight- 
binding model reproduces the ab initio band structure close Fermi energy with very good 
accuracy. 
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Table 5: Identification of the tight-binding 2 x 2 block structures and their orbital contributions 
at the r point with the band numbers and their corresponding DFT-HSE06 energies, given in 
Table 1^ In the last column, cases where the bands n, , and nj have the same energy at the 
T-point are denoted by nj-nj. 


block structure orbitals 

orbital composition 

(£+,£_) (eV) 

band numbers 

//£ : Oo — 03 


(-1.262,-7.571) 

(12VB, 6) 

He'. <Xi — 02 

P'x, Px 

(2.678,-4.105) 

(15-16,7-8) 

He (X2 — 01 

P‘y, Py 

(2.678,-4.105) 

(15-16,7-8) 

Ho'-Oti — 01 

dxz, Px, Px 

(2.457,-2.753) 

(13CB-14,10-11) 

Ho '.CCAr — 02 

dyz, Py, Py 

(2.457,-2.753) 

(13CB-14,10-11) 

Ho : 03 

PvPz 

-3.303 

9 


Although our tight-binding model contains many adjustable parameters, the optimization 
procedure presented in the Sec. [^imposes several implicit constraints. In practise, we find 
very difficult to obtain a parameter set that reproduces with high accuracy the position of 
the energy bands, their orbital compositions, and effective masses at the ^r,r, and Q points 
for both CB and VB. For this reason, in Table we present two parameter sets: one that 
reproduces most features of both VB and CB, but does not yield accurate masses for the VB, 
and the other that focuses on the VB. 

Table 6: Tight-binding model parameters obtained by optimization using a = 3.16A, 9b = 
0.710, and d = 2.406 A. The second column gives the best parameter set we obtain to fit 
both the valence (VB) and the conduction (CB) bands, while the third column focuses the 
optimization on the valence band. 


parameters 

CB-VB optimization (eV) 

VB optimization (eV) 

Do 

0.201 

0.191 

Di 

-1.563 

-1.599 

D 2 

-0.352 

0.081 

Dp 

-54.839 

-48.934 

Dz 

-39.275 

-37.981 

^pdic 

4.196 

4.115 

^pdo 

-9.880 

-8.963 

yppa 

12.734 

10.707 

^ppn 

-2.175 

-4.084 

^dda 

-1.153 

-1.154 

^ddn 

0.612 

0.964 

ydds 

0.086 

0.117 


Figure shows the tight-binding band structure calculated with the VB-CB optimized 
parameter set given in Table superposed with the DFT-HSE06 result. We find a very good 
agreement for the conductance and valence bands energies. The accuracy of the tight-binding 
results becomes increasingly poorer for band energies further away from the gap region, which 
is expected given that they were attributed a small weight in the fitting procedure. 
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Figure 5: Comparison between the band structures obtained with the DFT-HSE06 (blue) and 
with the optimized tight-binding model using the parameters from the CB-VB optimization 
(red) near the gap region. 


For completeness, in Table |7] we also include a comparison of the main orbital 
composition obtained from DFT-HSE06 and the tight-binding result. The orbital 
compositions at the high symmetry points r,K, and Q are not equal to those obtained with the 
DET-HSE06, but they show the correct leading d and p orbitals for both CB and VB. We point 
out that this is not the case in the tight-binding parametrization of Ref. Il25l . where several 
bands near the main energy gap have incorrect compositions. In particular, at the K point, 
the correct composition of the VB appears at a high-energy band, far from the gap. In our 
parametrization, out of the 18 points used in the optimization where analytical expressions 
where employed, only four yield incorrect compositions and they are located away from the 
main gap, at low energies. 

As shown in Table the effective masses are also reasonably well described by the 
CB-VB parametrization for all three special k points, except for the hole effective mass at 
the F point. We try to circumvent this limitation by performing another optimization (named 
VB) with a heavier weight given to the k values near symmetry points at the VB band. The 
resulting parametrization describes much more accurately the VB alone, imposing only few 
distortions on the CB, as Fig. [^reveals. This procedure yields the VB optimization parameters 
given in Tableand the orbital compositions and effective masses presented in Tables]^ and 
1^ respectively. Most parameter values are close to those of the global optimization (Table [^, 
while a few differ by more than 25%. The orbital compositions are nearly identical to those 
obtained in the global optimization. The most striking change is in the hole band mass the F 
point, which become essentically identical to the DFT-HSE06 value. 














































A tight binding model for M0S2 monolayers 


14 


Table 7: Top: Contribution from each orbital at the T point using the CB-VB optimization. 
Bottom: same but at the K point. Omitted orbitals have negligible or null contribution, t-b 
stands for tight-binding model. 


band number 







(T point) 

Py 

Pz 

Px 

dy^ 

d^2 

dxz 

DFT-HSE06 12 VB 
DET-HSE06 13 CB 

0.302 

0.141 

0.045 

0.327 

0.596 

0.049 

t-b 12 VB 

t-b 13 CB 

0.11 

n 

1 

0 

6.6 10-® 

0.889 

0.985 

5.4-10-® 


band number 
(K point) 

Py 

Px 

^xy 


dxi 

DET-HSE06 12 VB 

0.065 

0.065 

0.345 


0.345 

DET-HSE06 13 CB 

0.034 

0.034 


0.753 


t-b 12 VB 

2.7-10-4 

2.7-10-4 

0.499 


0.499 

t-b 13 CB 

00 

0 

1 

8.9-10-3 


0.982 



Table 8: Effective masses (in units of the free electron mass) at the T and K points resulting 
from the CB-VB and the VB optimization. 



HSE06 

CB-VB 

VB 

r point me 
r point m/, 

0.76 

-2.47 

0.35 

-0.62 

-2.47 

K point me 
K point;«/, 

0.42 

-0.47 

0.58 

-0.61 

-0.62 

Q point mg 

0.59 

0.59 



6. Simplified model 

The tight-binding model we have developed provides an accurate description of the main 
features of the CB and VB at the expense of involving a relatively large number of orbitals 
and fitting parameters. It has already been shown by Liu and coworkers ini that using 
just three orbitals for the Mo atom and including only the hopping amplitudes between Mo 
atoms in plane up to first neighbours is sufficient to open a band gap. With this in mind, 
we explored whether it is possible to neglect some hopping amplitudes in our tight-binding 
model and still obtain a reasonable description of the electronic structure near the band gap 
region. Keeping only the hopping amplitudes between Mo atoms turns out to be insufficient, 
as it preserves a large amount of degeneracy in the bands. Adding the hopping amplitudes 
between Mo and neighboring S atoms, without including the hopping amplitudes between S 
atoms, yields reasonable results. On the other hand, keeping exclusively the Mo-S hoppings 
does not yield a band gap. In matrix format, this simplified tight-binding model yields the 
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Figure 6: Comparison between the band structures obtained with the DFT-HSE06 (blue) and 
with the optimized tight-binding model using the parameters from the VB optimization (red) 
near the gap region. The VB optimization focuses on reproducing accurately the valence 
band. 

Table 9; Top; Contribution from each orbital at the F point using the parameters from the 
VB optimization. Bottom; Same but at the K point. Omitted orbitals have negligible or null 
contribution, t-b stands for tight-binding model. 


band number 



(F point) py 

Pz 

Px dy^ dyi dxz 

DFT-HSE06 12 VB 

0.141 

0.596 

t-b 12 VB 

1 .2-10-2 

0.988 


band number 
(K point) 

Py 

Px 

dxy d ^2 

dx2 

DET-HSE06 12 VB 

0.065 

0.065 

0.345 

0.345 

t-b 12 VB 

6.4-10-4 

6.4-10-4 

0.499 

0.499 


eigenvalue/eigenvector problem 

/ h^° + V r \ / a \ / a \ 

(py 0 T =£ T . ( 18 ) 

V {py 0 p J\p) \p J 

The number of fitting parameters is reduced to from 12 to 10. Symmetries can be 
fully exploited to break the diagonalization problem into smaller ones, as done previously. 
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After optimization against the DFT-HSE06 band structure, we obtain the values for the fitting 
parameters listed in Table 10 The resulting band structure is shown in Fig. [^superposed 
with the DFT-HSE06 band structure. We note that we were able to reproduce quite well the 
the entire VB, while the CB is well reproduced just around the K point, missing the correct 
behaviour aroung the Q and T points. Therefore, the simplihed model is somewhat limited 
in its applicability. It is suitable for the hole-doped region when the Fermi energy is brought 
to the top of the VB. It also provides a good description of the system when there is weak 
electron doping. 


Table 10; Parameters for the simplihed tight-binding model. 


parameter 

value (eV) 

Do 

-11.683 

Di 

-208.435 

Di 

-75.942 

Dp 

-23.761 

Dz 

-35.968 

Vpdp 

1.318 

Vpds 

-56.738 

Vdds 

-2.652 

Vddp 

1.750 

Vddd 

1.482 



Figure 7: Comparison between the DFT-HSE06 band structure (blue points) and the best ht 
to the simplihed tight-binding model (red points). 
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Due to the broken lattice inversion symmetry, strong spin-orbit interactions split the spin- 
degenerate valence bands in M 0 S 2 monolayers as well as in other group VI dichalcogenides. 
The spin-orbit coupling in this case is due to the Dresselhaus mechanism. Interestingly, the 
spin splitting in inequivalent valleys must be opposite, as imposed by time-reversal symmetry. 
As mentioned in Sec. [T] these features open interesting possibilities for the control of spins 
and valleys in these 2D materials IS a mill. 

Let us focus on the large spin-splitting at the K point of the VB. Its origin is qualitatively 
well understood: The valence band states are mostly made of dxy and dy 2 _y 2 orbitals with 
I = 2 and m = ±2. Therefore, the component of the SOC naturally gives a valley- 

dependent splitting of the bands. In contrast, the dominant contribution of the CB lowest 
energy state comes from the d^i orbital with I = 2 and m = 0, which cancels the spin-orbit 
splitting. These arguments agree with the quantitative analysis presented in Ref. IBl 

A complete tight-binding model that accounts for the effect of SOC over the entire 
Brillouin zone, including explicitly the p-orbitals of the chalcogen atoms, and taking into 
account the correct orbital composition of the main bands, is lacking. 

In this Section, we present an extension of our tight-binding model that includes the 
effect of an atomic spin-orbit coupling on all the atoms. For that purpose, we follow the 
formulation presented in Ref. Our starting point is the 11-band tight-binding spinless 

model derived earlier, with the Hamiltonian expressed in the appropriate symmetrized form, 
namely, where the block Hamiltonians He and Hq appear explicitly. The spin-orbit coupling 
term is inserted in the Hamiltonian by means of a pure intra-atomic spin-orbit interaction 
acting on all the atoms, explicitly given by 

^so = E^L,.S„, (19) 

where is the intrinsic effective SOC constant for an a atom (Mo o S), is the atomic 
orbital angular momentum operator, and Sa is the electronic spin operator. Hence, 


/ He 



( 

< 

<0 



Ho 

4 _ 




V 

He 




<0 

<E 


Ho y 


\ 

"ilr 



™oo 


The matrix elements of M are straightforward to obtain and depend on the SOC parameters 
Amo and As. The explicit form of the M matrices can be found in Ref. Il26ll . We note that in 
Eq. ( |20j i both diagonal and off-diagonal (spin-flip) terms are taken into account. However, an 
analysis in Ref. Il2^ indicates that spin-flip terms have a negligible contribution and could be 
dropped. 

We use DFT-HSE06 to estimate the splittings due to spin-orbit coupling and obtain 
= 202 meV at the K point of the VB. This value is higher than the experimental one ll47l . 
A^p = 145 ±4 meV. This is a known limitation of the HSE06 functional. The strength of 
this functional relies on its accuracy to predict the band gap of numerous materials, including 
M 0 S 2 , where traditional DFT calculations (EDA or GGA) give significantly understimated 
results. Using the SOC values Amo = 86 meV and As = 0.52 meV, which were obtained from 
a tight-binding parameter ht to maximally localized Wannier orbitals and to DFT calculations 
144), we hnd =173 meV. A better result is obtained by adopting the SOC parameters of 
Cappelluti et al. 12^ . namely, Amo = 75 meV and As = 0.52 meV. Inserting these values into 
our tight-binding formulation results instead in = 151 meV, which is in good agreement 
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with the experimental value. Thus we present in Fig. |^our results for the spin-resolved band 
structure based on this choice of SOC parameters. 



Figure 8: Comparison between the DFT-HSE06 spin-resolved band structure (blue points) 
and the best-fit tight-binding model (red points). The spin splitting is due to inclusion of 
spin-orbit coupling. 


8. Conclusions 

In conclusion, in this paper we provide a suitable and straightforward tight-binding model 
for a monolayer dichalcogenides, focusing our attention on M 0 S 2 . We show that this model 
reproduces rather well the structure near the main energy gap provided by an accurate DFT 
band structure calculation based on the HSE06 functional. It also reproduces the correct 
orbital composition of the bands. A fundamental ingredient in obtaining this result is the use 
of an optimization process that makes use of analytical expressions of the energy bands at 
symmetry points. In the constructing of our model we exploited the decoupling that exists 
between even and odd bands upon z inversion. Around the main gap, the bands are primarily 
even. Overall, the model yields 11 bands in the absence of spin-orbit coupling and involves 12 
fitting parameters. We provide two parametrizations for this case: one that is suitable for both 
conduction and valence bands (but less accurate for the valence band), and another that gives 
a very accurate description of important features of the valence band, such as the effective 
mass. When spin-orbit coupling is added, the number of fitting parameters jumps to 14. Our 
choice of parameters in this case yields a spin splitting of the valence band in good agreement 
with experimental values. 

We also investigate the possibility of turning off some hopping amplitudes in our model 
to reduce the number of parameters to 10 in the absence of spin-orbit coupling. The simplified 
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model is suitable for describing the hole-doped region or when one is only interested in the 
region around the K point. 

The present work provides a sound starting point for any further investigation of elec¬ 
tronic transport properties of single-layer semiconductor transition-metal dichalcogenides, or 
any other investigation that relies heavily on an accurate energy level positioning and wave 
function composition. 
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Appendix A. Tight-binding model and Slater-Koster parameters 


The tight-binding model is dehned by the following second-quantized Hamiltonian (spin 
indices have been omitted): 


r, 


'LL^V ifivY piv+{ptv,cyy Pl 


+H-C. 

r; H,v 

„t ,r,b A h 


'L L L f'llvdl^Pyy +t'llvdl^p'j,y +n.C. 

r; rj=r,+Ri-R2M.v 


ILL 

r; Tj=Ti-R2tl:V 


,l.t tI" t I jb tI [) I fj 

Vv dj i^ Pj^y + t^y dj^^ Pjy + H.C. 


'LL Wv'(p',v) pf,v'+H-C. 

r; v,v' ^ 

'L L \'’'L'dLdj.^i'+ii-C. 




+ LLi“vV {pivY p'jy + (ptv) p).y 
{ij)vy L 


H.C. 


(A.l) 


where {ij) denotes a sum over pairs of nearest-neighbour cells. The operators id] 
annihilate (create) an electron on the Mo in the unit cell i in the orbital p. Similarly, the 
operators [(pfy)^] and pi ^ [{pi y)^] annihilate (create) electrons at the bottom b and top 
t S sites of the unit cell i, respectively. We assume that the top and bottom S layers are 
symmetric (z inversion symmetry). 

We use the basis set defined in Sec. to express the on-site energies and the hopping 
integrals. The on-site energies are given by 

= (A.2) 


and 


gy = (r; + 5i±;pv|Jf'|r; + 5i±;py). 


(A.3) 
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The hopping matrix elements between Mo and S orbitals are 

= (r,;£/^|Jf’|r, + 5i+;pv), (A.4) 

f*v = (r,;t/^|^|r, + 5i_;pv), (A.5) 

= (r,;^^^|.^|r, +Ri -R 2 + 5i+;pv), (A. 6 ) 

f;;t = (r,; I Jf’lr; + Ri - R 2 + 5 1 _;pv), (A. 7 ) 

t’’[ = {ri-df,\Jif\ri-R2 + Si+-py), (A. 8 ) 

f^’t = (r,;fl'^i|^|r,-R 2 + 5i-;pv). (A.9) 

The hopping matrix elements between top and bottom S orbitals read 

Svv' = {ri + Si+-,p^\Jif\ri + Si--,py>), (A.IO) 

while the nearest-neighbor Mo-Mo hopping integrals are 

''mm' = {ri;d^\d^\ri + Rp,d^,), (ATI) 

''mm' = {ri;d^\^\ri + R2-,d^>), (A. 12) 

= {rr,d^\Jif\ri + R2-Rud^,), (A.13) 

'’mm' = {ri-,d^\J^\rt-Rv,d^,), (A. 14) 

''mm' = {ri-,d^\-^\ri-R2-,d^,), (A. 15) 

''mm' = {rr,d^\-d^\ri-R2 + Rud^,), (A. 16) 

(A. 17) 


and the S-S next-nearest-neighbor hopping matrix elements read 

“vv' = (ri + ^iiiPvl^lT+Ri+5i±;pv')> (A.18) 

“vv' = (T'+ 5i±;pv|'^|ri+R 2 + 5i±;pv')i (A.19) 

<^ = (r,- + 5i±;pvl^|r,+R2-Ri+5i±;pv'), (A.20) 

«vv' = (ri + 5i±;pv|-^|r/-Ri+ 5i±;pv'): (A.21) 

= (r/ + 5i±;pv|T^|r,'-R2-f 5i±;pv'), (A.22) 

“vv' = (*’! + 5i±;pv|>^|ri~R2 + Ri+ 5i±;pv')- (A.23) 


Notice that contains Mo-S and S-S nearest-neighbour (same unit cell) hopping 
amplitudes and Mo-Mo and S-S next-to-nearest-neighbour hopping amplitudes (adjacent 
cells). For the latter, each Mo and each S has six next-to-nearest neighbours. For the former, 
each Mo has six S nearest neighbours. Overall, there is a total of 25 hopping amplitudes 
within the unit cell and between the unit cell and the adjacent cells. The hopping amplitudes 


are indicated in Fig. A1 


r.MO 


The following associations are made for the on-site energies of Mo atoms: Cq ° = Aq, 
= e 


Mo _ 


= A 2 , £3 ° = £4 ° = Ai. For on-site energies of the S atoms we define £j = £3 = Ap 
and £3 = Aj;. 

The hopping amplitudes can be written in terms of SK integrals. We can also incorporate 
the X and z reflection symmetries, when applicable, to reduce the number of terms. Here we 
provide expressions for the more relevant hopping amplitudes in terms of seven SK integrals. 
This allows us to substantially reduce the number of fitting parameters of the model. Below 
we list amplitudes that are not identically zero. 
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Figure Al: Scheme of the hopping amplitudes. Solid black circles represent Mo atoms, while 
empty circles represent the S atoms at the top and bottom layers. 


Mo-S (Here we present the expressions for the t‘ hopping matrix elements. The t^ ones 
follow similar expressions, with 0* o —0*.) 


tQ2 = cos Oh ^sin^ ~ ^^hVpdTi, 

fo3 = sin 0 fo ^sin^ ^ ^ dtVpd,r, 

tn = — J- cos^ ObVpda - cos Ot sin^ ObVpdn, 


.2 .7 A 


f'lj = - sin Ob COS^ 9b ( —Vpda - Vpdn 1 , 


(A.24) 


(A.25) 


(A.26) 


(A.27) 


f21 = COS dbVpdx, 
f3i =&mebVpdn, 

Al = V^cos^0fosin0iyprf(j + sin0fo(l -2cos^0i)yprf;c, 

Al = sin^ 0focos 9bVpda+co& 0^(1 - 2sin^ 9b)VpdK, 

it ( 1 1 ? A ^ 1 

ffl’i = - — cos 0fo ( sin Ob--cos 0^, 1 y^^^j + - cos 0* sin OtVpdx, 


= 0*--cos 0fojyprf(j + —cos0fosin OhVpdx, 

43 = sin Ob ^sin^ ^pd(r + sin Bb cos^ ObVpdn, 


ti\ = - g COS^ 9bVpda - — cos 0fo 1^1 - - cos Ob j VpdK, 
^ ^bVpda + 7os 0* ^1 + ^ COS^ Ob^ Vpdn, 


(A.28) 

(A.29) 

(A.30) 

(A.31) 

(A.32) 

(A.33) 

(A.34) 

(A.35) 

(A.36) 
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t[i = ^ sin Ob cos^ Oh ^Vpda - Vpd^, 




'21 “ 


^22 = -^ QbVpda - ^ COS 0* I 1 - - COS^ ) Vpdjt, 


3 V 3 


cos'* ObVpda - - COS 0* 1 - - COS 0fo Vpdn 


V3 


1 


i,t V3 . 
f23 = ^sin 


j,t 3v^ 2 

hi = - 


Ob cos^ Ob [s/ZVpda - 2Vpdn 


COS Ob sin 9bVpda + sin 0^ 1 - - cos 0^ Vpdn, 


11 . 
%2 = 


db COtp- 9b {s/ZVpda - 2Vpdn 


3 "v/s 

f 33 = - - sin^ 9b cos 0foyprfcj - — cos 0^ (1 - 2 sin^ 9b) V, 


pdK') 


J,t . 

hi - 4 


sin 9b cos^ 0 fo {V^Vpda - 2 Vpdn 


^2 = ^ sin 9bVpda + sin 0fo ( 1 - cos^ 9b ) Vpdx 


V3 


■\/3 1 

f 43 = - ^ Sin^ 0* cos 9bVpda - - cos 0fo (1 - 2sin^ 0*) Vpdn- 


Mo-Mo (v'^, and s^^ can be obtained from and by symmetry.) 


22 

(A.37) 

(A.38) 

(A.39) 

(A.40) 

(A.41) 

(A.42) 

(A.43) 

(A.44) 

(A.45) 

(A.46) 


''00 — ^^dda + -^VddS, 

E V3 V3 
''01 —- -^ydda + —^VddS, 

E 3 1 

''11 — -^ydda + -VddS 


^22 — ^ddn, 

E T7 

^33 “ ^ddn^ 
1^44 = ^dd8 1 


'' 02 ^ — —Qydda + oVddS, 


(A.47) 

(A.48) 

(A.49) 

(A.50) 

(A.51) 

(A.52) 

(A.53) 


i,NE . 


..NE 


3v^ 


^/3, 


^/3, 


^ddoA- ^ Vddn ^^^ddS, 


V3, 


V3, 


''34 — ^Vddn - ^VddS- 


(A.54) 

(A.55) 
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• S-S («", m*", and can be obtained from and by symmetry.) 

£ ^ , 

^11 “ 'ppo, 

“22 ~ ^ppTll 

M33 = Vppji, 

U^2 = O^ppo ~ ^PP^) ■ 

■ 511 = ^ppn^ 

S22 — ^ppn^ 

S33 = Vppa- 

A general yet compact expression for all amplitudes is given by the matrices 


t = 


t^ = 


( 0 

0 

4i 

4i 

V 0 

( 0 

0 

t* 


ft 

'02 

^12 

0 

0 

ft 

M 2 


‘03 

t‘ 

'13 

0 

0 

t‘ 

'43 




-t 


31 


_ 


p = 


tpb ^ 


V 0 
/ 

i,, 

Ml 
/.r 
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/ 24 

'01 
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/.f 
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/ _ /4 
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\ -f'’' 

l,t 
01 
l,t 
'll 
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'02 
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13 
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‘02 
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hi 

j,t 

hi 

It 

'32 

j,t 
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J^t 
M 3 
l,t 

'33 

J^t 

'43 


/,r 

'02 

'12 

/,r 

'22 


_ \ 

‘03 \ 


-f 

-f 


‘03 
_ /.r 
M 3 
_ /4 

‘23 

l,t 

hi 

It 

hi 


( 


hi 

hi 

i,t 

i,t 

'12 

hi 

J,i 

ftL 

hi 

hi 

i,t 

i,t 

hi 

hi 

i,t 

i,t 

'42 

M 3 

i,t 

Lt 

hi 

hi 

i,t 

J.t 

hi 

-hi 


_ i,t 
“'22 
l,t 

hi 

_ j,t 

“M 2 


Lt 

hi 

Lt 

hi 

i,t 

'43 


(A.56) 

(A.57) 

(A.58) 

(A.59) 

(A.60) 

(A.61) 

(A.62) 


(A.63) 


(A.64) 


(A.65) 


(A.66) 


(A.67) 


(A.68) 
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(A.69) 


(A.70) 




vA'^ = 




4^00 

2^01 

NE 

^02 


0 

0 

\ 


_ lyfi 
2^01 

IE , 3 E 

4^11 + 4422 

NE 

^12 


0 

0 



4 o 2 

NE 

4 i 2 

4 V 22 + 4^11 

0 

0 


, (A.71) 

0 

0 

0 


4vf3 + |vf4 

NE 

V34 



0 

0 

0 


< 

iy£ + 3 £ / 

44’44+ 4V33 / 


^00 

_iy£ 

2^01 

_yA'£ 

^02 


0 

0 




^vfi + |vf2 

_yA'£ 

^12 


0 

0 



-vS'# 


1 £ ,3, 
4^22+ 4 

4^11 

0 

0 


, (A.72) 

0 

0 

0 


iy£ _L 3 £ 

4 V33 -h 4 V44 

-< 



0 

0 

0 


_y'V£ 

434 


4’f3 ) 



wfj 0 

0 

\ 



0 M22 

0 

. 


(A.73) 

0 0 

1*33 

) 



44*11 + 4 

**f2 




**12 


1, .£ ,3 , .£ 

41*22+ 4**11 

0 , 

(A.74) 

0 


0 

1*33 / 





<11 + 4M22 


(A.75) 


Appendix B. Tight-binding energy bands 

To find the energy bands we need to solve the the eigenvalue/eigenvector problem in the Bloch 
momentum representation, 

J^"|k)=£(k)|k), (B.l) 

where the Bloch vector is given by 

[43 

|k> = ^ (j3k.|j|ri + 5i_;p^) + Tk,,j|ri + 5i+;p^)) . (B.2) 

r, /j =0 11=1 

Resorting to the orthogonal bais, we can rewrite the eigenvalue/eigenvector problem as a 
system of linear coupled equations, 

rctrt _ cMol _ V r A' _L T* _L T* T. (■'■4 3 ,r. _L T* _L T* T. R. 1 


[^(k) -e“°] ttk,^ — (f[iv +4f[n/ +Z 2 ^d'ifv) 'fk,v + 

+ 2^(civ^^,+C2vJ^|j,+Cl2vJ^^^ Ctk.Ai', 
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and 


[£(k) - gS] j3u,v = E {t% + Z2t‘fA +Z 2 Z*/ff,) au,^ 

+ H [(*^1 “vv' +‘^2Mvv/ +Ci2Mv^) Pky +«VV' '^k.v'] ■ 

v' 

[£(k) - e^] Tk,v = ^ (f^v 

+ H [(‘'1 “vv' +‘22Mvv' +C12 Mv*) +^vv’Pky] j 
•»»/ 


(B.4) 


where zi 


Cl = 


and 


gik-R, _ g!a/:;c _ gik-R2 _ ^iak^jl^ifhakyjl 

= = COS(k• Ri) = cos (kxa), 

C 2 = =cos(k-R 2 ) =cos (j(xa/2. + ky\/3a/2^ , 

cj 2 = ^2 ^2 _ • Ri — k • R 2 ) = cos {kxa/2. — ky\/5a/2}j 


In matrix form, we have 
/i“°+y 

{rf h^+u 

/ . s T 


{t'’) s h^ + U 


a 

X 


a 

= E\ I |, 

(3 


(B.5) 

(B.6) 

(B.7) 

(B.8) 

(B.9) 


where 


h^° = 


/ Ao 0 0 0 0 \ 

' 0 A2 0 0 o' 

0 0 A 2 0 0 

0 0 0 Ai 0 

\ 0 0 0 0 Ai y 


. Ap 0 0 

/i^ = I 0 Ap 0 


0 0 A, 


r = f+zl{t^^‘+zi 






and 


+ zl{t‘^‘’ + zif 
V = 2 (ci + C2 + C12 V 


,NW\ 


(B.IO) 

(B.ll) 

(B.12) 

(B.13) 

(B.14) 


t/ = 2(ciM^+C2M'^‘^ + Cl2M^'^). (B.15) 

This formulation suffices for a numerical evaluation of the bands. However, in order 
to obtain analytical expression for the gap and other features of the band structure at the 
symmetry points, it is necessary to reduce the size of the matrices to be diagonalized. This 
can be done by exploring underlying symmetries in the equations. 
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Let us define the symmetric and anti-symmetric components 


^,v = (Tk.v + Pk,v) 


(B.16) 


•/•k.v — ("^kjV ~ ^k,v)- 

In terms of these components, the eigenproblem takes the form 

/ h^° + V \ / C 

h^ + U + s 0 ( 

V 0 h^ + U-s j \ (I 

where we introduced new hopping matrices 




0 =£ 0 


(B.17) 


(B.18) 


(B.19) 




(B.20) 


We can further simplify the eigenproblem by rearranging amplitudes in the eigenvector, 
going from 

V/^ = (ao,ai,a2,«3,«4,0i,02,03,^i,02,03)- (B.21) 


= (ao,ai,a2,01,02,03,«3,«4,01,'/>2,03)- 


(B.22) 


This amounts to ordering the basis states such that the first six components in the eigenvector 
are even (E) while the last five are odd (O) with respect to z inversion. As a result, the 
eigenproblem can be recast in the decoupled form 


where 


0 

0 


W = EW, 


^Mo.E "yE^O 

rfi.ot hS,E 


(B.23) 


(B.24) 


\ rO.£t j ■ 

We have introduced the following matrices: 

/ Ao 0 0 \ 

/jMo,E^ 0 A 2 0 +y^, 

VO OA 2 


/jMo ,0 ^ 


Ai 0 

0 Ai 


Ap 0 0 


0 Ap 0 + 

\ 0 0 A, / 


(B.25) 


(B.26) 


(B.27) 


(B.28) 
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and 


/ Ap 0 0 \ 

0 Ap 0 + 

Vo 0 aJ 



/ Voo Vbi V 02 \ 

V^=\ Voi Vii V12 , 

V V02 Vl2 V22 / 


( V 33 y34 V 

V ^34 V 44 ) ’ 


/ Un Ui2 0 \ 

Vo 0 t/ 33 y 

/ -ill 0 0 \ 

= 0 -*22 0 , 

Vo 0 *33 / 


/ ill 0 0 \ 

5 ^ = 0 *22 0 

Vo 0 -*33 J 


(B.29) 

(B.30) 

(B.31) 

(B.32) 

(B.33) 

(B.34) 

(B.35) 

(B.36) 

(B.37) 


Appendix C. Expansion around symmetry points 

Expanding the hopping matrix elements around symmetry points in the Brillouin zone allows 
to obtain analytical expressions for bands energies and orbital composition. 


Appendix C.l. T point 

At the r point, kx = ky = 0, resulting in zi = Z 2 = 1 and ci = C 2 = C 12 = 1- Then, 

/ 

0 


rE,0 


= V2\ 


0 

0 

ti, + 2f. 


Lt 


ft _j_ J-t 
(l2-|-,(i2 

0 


f03 +2fQ3 

0 


0 


7^*’^ = V 2 


?3i H-2t 


= 2 


34 

0 

0 


142 + 2 ? 


42 


5 (v^ 
2 ril 


'' 22 ! 


i(vfi+- 


22 ) 


(C.l) 


(C.2) 

(C.3) 
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and 


‘ 02 ' 


‘02 


Vf = 2 


Ur = 2 


Sr = 


Sr = 


'13' 


(1 

(vf3 + vf4) 


\ 

0 



K 1 +M 22 ) 



0 


V 

0 


-■Si 

1 0 

0 

0 

—S22 

0 

0 

0 

S 33 

Sll 

0 0 


0 

S22 0 


0 

0 -S 33 y 

l.t 

13 — 

^43 +2f43 = 

= 0 


vfj- 


2 V '^11 


'' 44 } 


^22) 


0 

0 

3 M 33 


(C.4) 


(C.5) 


(C. 6 ) 


(C.7) 


into three 2x2 diagonal blocks and into two 2x2 blocks and one 1x1 block. All blocks 
can then be diagonalized analytically. Explicitly, 


• 00 — 03 


• 01-02 


• 02-01 


• 0 : 3-01 


• 04 — 02 


• 03 


/ Ao + 6 v^o 

y •\/2 +2fQ3) 

/ A 2 + 3 (vfi Avfj) 

72(4 + 2 ?;^) 

/ A 2 + 3 (vfi Avfj) 

v^( 4 i+ 24 ') 

/ Ai +3 (v33+vf4) 
/ Ai +3 (v33+vf4) 

^2(4+24^) 

( Aj + 6 Mf 3 +S 33 ) . 


s/2 ^fQ3 +2fQ3^ 

Az + 6Mf3-i33 


(^12 + ^^12 
Ap + 3 (wfi + U 22 ) 

3/2 + 2 f 2 j 

Ap + 3 (wfi + Mfj) 

y/2 (^2,\ + 2 f 3 \ 
Ap + 3 (wfj +U 22 ) 






Ap + 3 (wfi + wfj) 


— S 22 


(C. 8 ) 

(C.9) 

(C.IO) 

(C.ll) 

(C.12) 

(C.13) 


The bands can be identibed by matching their composition to the results of the DFT- 
HSE06 calculations, (see Table [^. For instance, the valence band state at the E point is 
mostly composed by d ^2 and pj.-orbitals. Another band state (band number 6 , as debned in 
Table 1^ that has a similar orbital composition is lower in energy than the valence band state. 
Therefore, we associate the highest eigenvalue of (|C. 8 |l to the valence band energy, while the 
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lowest eigenvalue is put in correspondence with the band number 6 energy. As a result, we 
find 


^valence (r) — 


Ao + 6 Vaa + Aj; + 6 wfo — S33 


(Ao + 6vfo-Aj-6Mf3+S33)' 


+ 2 ^to3 + 2 f/’^ 


03 


(C.14) 


and 


£ 6 (r) = 


Aq + 6 Vqq + Aj + 6 m?3 — S33 


(Ao + 6 vgo - Aj - 6 Mf3 + S33) 


(^03 


+ 2 (f^3+2fo’3 


(C.15) 


Carrying out the same procedure for other blocks, we arrive at fully analytical expressions for 
all bands at the F point. 

Thus each block correspond to a majority orbital composition and each eigenvalue (two 
for each block 2 x 2 £_ and F+) is matched to a DFT value. 

Appendix C.2. K point 

At the K point, kx = 2;r/3fl and ky = —2%l\/3a, resulting in zi = Z 2 = —and 

Cl = C 2 = C 12 = -1/2. Then, 


rE,0 


■.0,E 


= V2\ 


= V2\ 


is/3t 


^31 


'01 

d - 
'02 

/,r 

'02 

l,t 

'11 

d - 
'12 

l.t 

^12 

l,t 

'21 

/ ?22 

t‘'‘ 

'31 

i 0/3/22 

J^t 

hi 

ff _ 

^42 

i,t 

'42 


0 

'13 '13 

i\/3t22 


iV3t 

^43 ~ ^43 


l,t 

33 

l.t 


(C.16) 


(C.17) 




0 

0 


0 

^fl+''f 2 

0 


0 

0 

4+''fi 


(C. 18 ) 


I 

II 

0^4 


vf3+vf4 

0 

V44 + vf3 ) ’ 


(C.19) 

b 

II 

3 

0 

/ My + U 22 
0 

0 

U 22 + My 

n. 

(C.20) 


0 


2 M 33 


cO _ 

Ot- — 


-ill 0 0 

0 —*22 0 

0 0 S 33 


(C.21) 


and 


cE _ 
Ok — 


ill 0 0 

0 S 22 0 

0 0 —*33 


(C.22) 


where we have imposed tQ^, = ^03 sn = S 22 - 
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The Hamiltonian matrix can be block diagonalized by a chiral transformation 

«! ± ia2 


(XR2,L2 = 


OtRlMl = 
Orx = 


V2 ’ 

as ± ia^ 
01 ± /02 


= 


V2 ’ 

01 ± i02 


(C.23) 

(C.24) 

(C.25) 


(C.26) 


The other variables, ao, 03 , and 03, remain the same. Thus, we have the new state vector 


V = (O!O,O:L2,O:«2,0L,0«,03,O!Ll5Ct«l,0L,0fi,03)- 


The result is 


H = 


(C.27) 


(C.28) 


where 


A = 


^o+Vqo 

0 

0 

0 

—i{Koi —K02) 

0 


0 

A2 + Vi^ 

0 

— —-^^21 +-^^12 + -^^ 22 ) 

0 
0 

0 

—K 21 +K 12 +K 22 ) 

0 

■Sll 


0 

0 

A2 + V22 
0 
0 

^13 — ^23 


(C.29) 


and 


B 


Al +^3*3 
0 
0 
0 

V -*(-^33--^ 43 ) 


-Un 

0 

0 




KKqi —K02) 
0 
0 
0 

Ap + U 22 + S22 
0 


0 
0 

K 12 —K 23 
0 
0 

+ 1 / 33-^33 / 


0 

Ai +V^ 

— A'41 — ■S'32 — -^42) 
0 
0 

0 


(C.30) 


(/Tai — iir 4 i — ^^32 — K42) 

Ap + t/ 11 —ill 
0 
0 

The following elements have been introduced; 


0 

0 

0 

Ap + U22 — ifii 


KK 23 -K 42 ) \ 
0 
0 
0 

AZ + U 32 + S 33 j 


Koi = y3foP 
Kn=V3t[’l, 
K 21 —(‘21 — t!'' 
■^31 = ^31 

Ku = y3f'’[, 


21 ’ 

l,t 

31’ 


Ko2=to2-to2^ 
K\2 =fl2-fl2’ 
K 22 = V^t22 , 
K22 = V3t: 


32’ 


IJ 


K42 — f 42 — f 42 , 


■^13 =tl3-hi’ 
K 23 = V3t!^i, 
K33 = v/ 3 f 33 , 
■K 43 = ^43 ~ ^43- 
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Using the SK decomposition, one finds that several combinations of these coefficients yield 
zero. These simplications have already been implemented in Eqs. ( C.30| l and (C.311. The H 
matrix breaks up into five 2x2 blocks and one 1x1 block, 


• OCq — Or: 


• an — Ol'. 


0tR2 — 03: 


• CCrI — <1>L' 


• CClI — 63 : 




^0 + ^00 i{KQi—Ko2) 
—i{Koi—Ko2) Ap + Ufj+iii 


A2+y/i 

—i{Kii —K21 +K12+K22)/V2 


A 2 + 1^22 

K 12 —K 23 


Ki3 —K 23 

Aj +1/33 i 33 


(C.31) 

*(■^11 —-^ 21 +-^ fl 2 +-^ f 22)/'\/2 \ 

^p + U 22 +Sll ) 

(C.33) 


(C. 


Ai+yo 

{Ki 2 — K 41 — K 22 — K 42 )/Vl 


{Ki 3 — K 41 — K 23 — K 42 )/V 2 \ 


An -f 1^22 1 




(C.34) 


Ai +y3*3 
-i{K22 - K42) 


Ap + t/j[ — ill 


i{K33-K43) 

JO- 


^z + t^33 +'S'33 


(C.35) 


(C.36) 


According to the DFT-HSE06 calculations, at the K point, the conductance band is 
mainly composed by d 2 ^ 2 _^i, Px, and py orbitals, while the valence band is mainly formed 
by d,2_y2, dxy, Px, and py orbitals. Therefore, the conductance band can be obtained from 
04) —Ol variables, while the valence band comes from the aL2 — Or combinations, resulting in 
the expressions 

"^0++Ap+{/^[ +iii 


E'conductance(^) — ‘ 


'(Ao+yo^o-Ap-t/fi-iii)2 


(i^oi —K() 2 Y 


(C.37) 


and 


7:valence(^) — 


A 2 +yn + Ap + f/^ +ill 


+ 


(A 2 + y/i - Ap - t/ f2-iii)2 ^ {Kn- K 21 - Kn - K 22 f 


(C.38) 


4 2 

Following the same procedure for other blocks, we find analytical expressions for nearly 
all energies at the K point. Thus, each block corresponds to a majority orbital composition 
and each eigenvalue (two for each block 2 x 2) is matched to a DFT-HSE06 value. 
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Appendix D. Alternative implementation of unsymmetrized band equations and 
comparison with the tight-binding model of Cappelluti et al. 


The essential difference between our construction of the 11-orbital and that of Cappelluti et al 
Il25ll2^ comes from their inclusion of two phase factors in the Bloch state equation, namely, 




4 

11=0 


ri\da 


M=1 




+Tk.^e'“"'+|r; + 5i+;p^)) . (D.l) 

Comparing this equation with Eq. we notice the extra phase factors in the amplitudes 

of the S atomic orbitals. While the phase factors have no impact on the eigenvalue secular 
equation, they do change the matrices containing hopping amplitudes between Mo and S 
atoms. For instance, our T' matrix of Eq. (B. 121 would change to 


r = 




+ Z\t 


r,t 


Jk-8 


'1 + 


(D.2) 


A second yet important difference between their work and ours is on the notation and 
organization of the hopping matrices. 

To facilitate a direct comparison between our model and that of Refs. Il25ll2^ . we begin 
by swapping the second and third rows and corresponding columns in Eq. ( |B.9| l, 

/ h^ + U {T‘y s \ / T \ / T \ 

r h^° + V T‘’ a =£ a . (D.3) 

V 5 h^ + U ] \ P ) \ P ] 


Next we introduce their auxiliary quantities ^ = kxa/2 and rj = s/^kya/2, which allows us to 
rewrite the coefficients in Eqs. ( [B^] ), ( [B^ , and (|^ as 

ci=cos(2^), (D.4) 

C 2 =cos(i*-ft]), (D.5) 


and 

C12 = COs(^ -T]). 

Also, zi = e^‘^, Z 2 = and 

g,k.5, ^ ^2in/3 
gi'k-52 ^ 
gi'k-53 ^ 


(D.6) 

(D.7) 

(D.8) 

(D.9) 


The correspondence between our block matrices and theirs is the following (phase factors 
set to zero in the appropriate hopping amplitudes); 

• h^ + U ^ Hpt^p, = Hpb^pb with 


Hpt,pt — 



• 5 O Hpt^ph, with 


Hpt,pb — 


H, 

H, 


xy 


yy 

0 


0 

0 


H-, 


zz,z 


V, 


ppn 

0 

0 


V, 


PPK 

0 Vr 


ppa 


(D.IO) 


(D.ll) 
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• h^° + y <->■ with 


Hd,d = 


T‘ o Hd^pt, with 


^d.pt — 


/ 

^z^z^ 

Hziyi 

^z^xy 


u* 

^x^x'^ 

Hx^xy 


H*2 

z xy 

H*2 

x'^xy 

^xy^xy 


0 

0 

0 

V 

0 

0 

0 

/ 

Hz2x 

^zh 

Hz^z 


^x^x 

Hxly 

f^xh 


Hxy^x 

^xy,y 

Hxy,z 


Hxz,z 

Hxz,y 

Hxz,z 


Hyz,x 

^yz,y 

^yz„z 




H- 


xz,xz 


H, 


xz,yz 


xz,yz 




yz,yz 


(D.12) 


(D.13) 


<->■ Hd^pb, with 


f^d.pb — 


Hz^x 

Hzh 

-Hzh 

^x^x 

H^ly 

-Hx2z 

^xy,x 

^xy,y 

~Hxy,z 

-Hxz,z 

^xz,y 

Hxz,z 

-Hyz,x 

~^yz,y 

Hyz,z 


V 


(D.14) 


All matrix elements are identical to those of Ref. 
Hd,pb, which are explicitly dehned below: 


/ 

, except for the matrices Hd^pt and 


= - 2 V 3 E 1 sm{^)di (D.15) 

H^2y =2 EiC2 (D.16) 

H. 2 ^ =£2 Cl (D.17) 

H^2, = - 2 V 3 Qcs-Ca) sin(^)r/i (D.18) 

Hxy,y=H^2, (D.19) 

Hxz,z = -2 sin(^)ifi (D.20) 

Hyz,x = H,z,y (D.21) 

Hy,^,=2E^C2, (D.22) 

H^2y = -2£3C3+2i£5Cos(^)<ii (D.23) 

H^2^ =2EaC2 (D.24) 

Hxypc = - ^CsCs + 6!£'3 cos(^)£/i (D.25) 

Hxy,z = —2y/3E4 sin(^)£/i (D.26) 

Hx,,y = - 2^3 Q £6 - C 7 ) sin(^) di (D.27) 

Hxz,x = - ^£ 6 C 3 +6^£7Cos(^*)^fl (D.28) 

Hy,^^y = - 2 £ 7 C 3 + 2iEb cos(^ )di. (D.29) 
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In all these equations the quantities h, I 2 and I 3 , as well Et, i = 1,..., 8 , follow the definitions 


of Ref. II 25 I : for instance. 

/i = cos(2^)+2 cos(^) cos(t 7 ) (D.30) 

/2 = cos(2^) — cos(^) cos(t 7 ) (D-31) 

/3 = 2cos(2i*)+cos(^) cos(tj). (D.32) 

The coefficients Ci, C 2 , C 3 , and di become more compact without the inclusion of phases, 
namely, 

Cl = — 1 —2 cos(^) cos(t 7 ) +2/cos(^) sin(T 7 ) (D.33) 

C 2 = — 1 +cos(^) cos(t 7 ) — icos(^) sin(77) (D.34) 

C 3 = cos(^) cos(t 7 ) — /cos(^) sin(77) +2 (D.35) 

t/i = !cos(t 7 )+ sin(T 7 ). (D.36) 
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